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PROPERTIES OF QUASI-RELAXED STELLAR SYSTEMS IN AN EXTERNAL TIDAL FIELD 
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ABSTRACT 

In a previous paper, we have constructed a family of self-consistent triaxial models of quasi-relaxed stellar 
systems, shaped by the tidal field of the hosting galaxy, as an extension of the well-known spherical King 
models. For a given tidal field, the models are characterized by two physical scales (such as total mass and cen- 
tral velocity dispersion) and two dimensionless parameters (the concentration parameter and the tidal strength). 
The most significant departure from spherical symmetry occurs when the truncation radius of the corresponding 
spherical King model is of the order of the tidal radius, which, for a given tidal strength, is set by the maximum 
concentration value admitted. For such maximally extended (or "critical") models the outer boundary has a 
generally triaxial shape, given by the zero-velocity surface of the relevant Jacobi integral, which is basically 
independent of the concentration parameter. In turn, the external tidal field can give rise to significant global 
departures from spherical symmetry (as measured, for example, by the quadrupole of the mass distribution of 
the stellar system) only for low-concentration models, for which the allowed maximal value of the tidal strength 
can be relatively high. In this paper we describe in systematic detail the intrinsic and the projected structure and 
kinematics of the models, covering the entire parameter space, from the case of sub-critical (characterized by 
"underfilling" of the relevant Roche volume) to that of critical models. The intrinsic properties can be a useful 
starting point for numerical simulations and other investigations that require initialization of a stellar system in 
dynamical equilibrium. The projected properties are a key step in the direction of a comparison with observed 
globular clusters and other candidate stellar systems. 
Subject headings: globular clusters: general — stellar dynamics 



1. INTRODUCTION 

It is commonly thought that globular clusters can be de- 
scribed as stellar systems of finite size, with a truncation in 
their density distribution determined by the tidal field of the 
hosting galaxy. Given the fact that many globular clusters un- 
dergo significant relaxation, the simplest analytical models 
( lKindll966h that provide a successful description of these 
stellar systems are defined by a quasi-Maxwellian distribu- 
tion function fxiE), where E is the single-star energy, with a 
truncation prescription continuous in phase space. 

Most of the interesting physical mechanisms that under- 
lie the dynamical evolution of these stellar systems (such 
as evaporation and core collapse; e.g., see ISpitzeii Il987t : 
iHeggie & Hutll2003l) depend on such truncation and are fre- 
quently studied in the context of spherical models for which 
the action of tides is implemented by means of the existence 
of a suitable truncation radius (which, in principle, for a glob- 
ular cluster of given mass can be determined a priori if we 
know the cluster orbit and the mass distribution of the host- 
ing galaxy) supplemented by a recipe for the escape of stars. 
Therefore, evolutionary models that rely on the assumption 
of spherical symmetry, such as Monte Carlo models ( e.g., for 
a description of two of the codes currently used, see iGierszl 
119981: iJoshi et allllOOOl) and Fokker-Planck models (e.g., for 
an exam ple of, respectively, the isotropic and anisotropic 
case, see IChemoff & Weinberg 1990; Takahashi et al. 1997) 
are necessarily based on an approximate treatment of the 
tidal field. In this context one interesting exception to the 
adoption of spherical sym metry is that of the axisymmetric 
Fokker-Planck models by lEinsel & Spurzenj (fT999l) . which 
admit flattening induced by internal rotation. 

Yet, if tides are indeed responsible for the truncation, they 
should also induce significant deviations from spherical sym- 
metry: in the simplest case of a cluster in circular orbit 



about the center of the hosting galaxy, the associated (steady) 
tidal field is nonspherical and determines an elongation of 
the mass distribution in the direction of the center of mass 
of the hosting galaxy accompanied by a compressio n in the 
direct i on perpendicular to the orbit plane (e.g., see ISpitze^ 
11981: iHeggie & Hutll2003h . Direct N-body simulations, in 
which an external tidal field can be taken into account ex- 
plicitly, provide a unique tool for the study of the evolution 
of a tidally perturbed cluster, especially when non-circular 
orbits are considered so that tidal effe cts are time-dependent 
(e.g., see iBaumgardt & Makinoll2003h ; in particular, this ap- 
proach has led to detailed models of the rich morphology of 
tid al tails, i.e. strea ms of stars escaped from the cluster (e.g., 
seelLeeetaT 2006). 

Indeed, deviations from spherical symmetry are observed in 
globular clusters (e.g.. see fOeyer et al]ll983l : rWhite & Shawll 
11983), but, they are often ascribed to other physical ingredi- 
ents, such as internal rotation. In fact, it is recognized that 
the issue of what determines the o bserved shapes of globular 
clusters remains unc lear (e . g.. see |King|ll96lHFall & FrenkI 
119851; lHan&Rvdenlll994l; Ivan den Berghl 120081; and refer- 
ences therein). 

One might argue that real globular clusters are likely to be 
not fully relaxed, may possess some rotation and experience 
time-dependent tides so that analytical refinements beyond 
the spherical one-component King models would not compete 
with the currently available numerical simulation tools (see 
§ 5.1 for additional comments) that allow us to include these 
and a great variety of other detailed effects that are relevant 
for the quasi-equilibrium configurations. However, physically 
simple analytical models, accompanied by the study of more 
realistic numerical simulations, serve as a useful tool to inter- 
pret real data and to provide insights into dynamical mech- 
anisms, even though we know that real objects certainly in- 
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elude features that go well beyond such simple physical mod- 
els. In this specific case, we argue that the triaxial geometry 
is a natural attribute of the physical picture of tidal truncation, 
which has already proved to provide useful guidance in the 
study of globular clusters. 

As demonstrated in a previous paper (iBertin & Varrill2008t 
hereafter Paper I), a complete characterization of the physi- 
cally simple description of globular clusters as quasi-relaxed 
tidally- truncated stellar systems, in terms of fully self- 
consistent triaxial models, can be provided. In order to better 
address the role of tides in determining the observed structure 
of globular clusters, in Paper I we have thus constructed self- 
consistent non-spherical equilibrium models of quasi-relaxed 
stellar systems, obtained from the spherical case by includ- 
ing in their distribution function the effects associated with 
the presence of an external tidal field explicitly. We recall 
that our models consider the stellar system in circular or- 
bit within the hosting galaxy, for simplicity assumed to be 
spherically-symmetric. Therefore, in the corotating frame 
of reference, the tidal field experienced by the system is 
static and the Jacobi integral H is available. In this physi- 
cal picture the typical dynamical time associated with the or- 
bits inside the cluster is assumed to be much smaller than 
the external orbital time. The procedure described in the 
previous paper starts by replacing the single-star energy E 
with the Jacobi integral in the relevant distribution function 
fK{H) = A[exp(-fl//)-exp(-fl//o)] if H < Hq, with Ho the 
cut-off constant, and fxiH) = otherwise. Thus the collision- 
less Boltzmann equation is satisfied. The construction of the 
self-consistent models then requires the solution of the asso- 
ciated Poisson-Laplace equation, that is of a second-order el- 
liptic partial differential equation in a free boundary problem, 
because the boundary of the configuration, which represents 
the separation between the Poisson and Laplace domains, can 
be determined only a posteriori. The idea of using the Ja- 
cobi integral for the construction of tidal triaxial models had 
been proposed also by lWeinbergI (119931) . A family of models 
based on the distribution function fK(H ), but constructed with 
a diffe rent method, was discussed by iHeggie & Ramamanil 
(fT995h : a comparison with our models will be given in § 3.2. 

Here we describe the properties of the resulting two- 
parameter family of physically justified triaxial models, con- 
structed with the method of matched asymptotic expansions 
(see § 4 and § A in Paper I), inspired by previous s tudies of the 
formally similar problem of rotating poly tropes ( Smith|[T975l; 
[l976). In particular, we illustrate the properties of the first 
and second-order solutions. 

The paper is organized as follows. In § 2 we present a thor- 
ough description of the relevant parameter space. The intrin- 
sic and projected density distributions are discussed in § 3, 
with special emphasis on the global and local quantities that 
can be used as diagnostics of deviations from spherical sym- 
metry. Intrinsic and projected kinematics are addressed in § 4. 
The concluding § 5 gives the summary of the paper with a dis- 
cussion of the results obtained and a comment on the complex 
physical phenomena that a large body of evolutionary models 
based on numerical investigations has shown to characterize 
the periphery of globular clusters. 

2. THE PARAMETER SPACE 

The triaxial tidal models are characterized by two physical 
scales (corresponding to the two free constants A and a in the 
distribution function fK{H)) and two dimensionless parame- 
ters. The latter parameters are best introduced by referring to 
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Fig. 1. — Parameter space for second-order models. The uppermost solid 
line represents the critical values of the tidal strength parameter for models 
in which the potential of the hosting galaxy is Keplerian (u = 3); thin solid 
lines represent contour levels of the extension parameter 5. The dashed line 
is the critical curve for models within a logaiithmic potential {u = 2). The 
dot-dashed line gives the critical condition for a potential with u = I (e.g., 
that of a Plummer sphere evaluated at Rq = h/\/2, with h the model scale 
radius). 

the formulation of the Poisson equation in terms of the dimen- 
sionless escape energy 

^(r) = aHo-[a<^>c(r) + eT(x,z)], (1) 

where fl$c is the dimensionless cluster mean-field poten- 
tial (to be determined self-consistently) and T{x,z) = 9{z^- 
vx^) /2 represents the tidal potential (with the numerical co- 
efficient V = 4-K^/f2^, where n and f2 are respectively the 
epicyclic and orbital frequency, depending on the potential 
of the hosting galaxy); the hat on the spatial coordinates 
denotes that they are measured in units of the scale length 

ro = [9/(4^Gpofl)]'/2. 

Then the first parameter, already available in the spherical 
case, is the concentration of the system and can be expressed 
as a dimensionless measure of the central depth of the poten- 
tial well: = V'(O). The second parameter, the tidal strength 
e in Eq. ([T]i, is defined as 



i.e. as the ratio of the square of the orbital frequency (of the 
revolution of the stellar system around the center of the host- 
ing galaxy) to the square of the dynamical frequency associ- 
ated with the central density po of the stellar system. Alter- 
natively, the effect of the tidal field can be measured by the 
extension parameter 

5 = ftr/h , (3) 

where r,r = r,r('i') is the truncation radius of the spherical 
King model characterized by the same value of ^ and i-j is 
the tidal (or Jacobi) radius, i.e. the distance from the origin 
(the center of the stellar system) of the two nearby Lagrangian 
points of the restricted three-body problem considered in our 
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simple physical picture. A given model will be labelled by 
the pair of values (5',e) or, equivalently, by the pair (^,(5). 
The dimensionless cut-off constant qHq can be expressed in a 
natural way as an asymptotic series with respect to the tidal 
parameter aHo = ao + aie+aze^ /2 + ... where the terms a/, as 
discussed in Paper I, depend only on ^. 

Much like the Hill surfaces for the standard restricted three- 
body problem, we now consider the family of zero-velocity 
surfaces defined by the condition 4'{r) = 0, which represents 
the boundary of our models. These surfaces can be open or 
closed, depending on the value of the cut-off constant a/Zo, 
which is determined by the selected values of the two dimen- 
sionless parameters that characterize the model. To be con- 
sistent with the hypothesis of stationarity, we only consider 
closed configurations. We call "critical models" those that 
are bounded by the critical zero-velocity surface (which is the 
outermost available closed surface). For each value of ^, the 
critical value of the tidal parameter can be found by (numeri- 
cally) solving the system 



d^iipix = rj- ,y = 0,z = 0; e„) = 
V'(i = rr,}' = 0,z = 0;e„) = , 



(4) 



where the unknowns are rj and €cr- The method of matched 
asymptotic expansions proposed in Paper I for the solution 
of the relevant Poisson-Laplace equation requires an expan- 
sion in spherical harmonics, therefore it can be easily recog- 
nized that the first condition of Eq. (|4|i is equivalent to the 
requirement of vanishing gradient, which identifies the saddle 
points of the critical surface. In the general case, the condition 
dxtpi^T, 0,0; e) = determines the value of rr for a given tidal 
strength e , therefore rj- = r7-(5', e). 

By using in the escape energy defined in Eq. ([U the zeroth- 
order expression for the cluster potential (fl$^"')*"'(r) = Ao/r 
and for the cut-off constant', the system in Eq. dUi becomes 



■ 



Ao 

'T 

Oio- — + -ecr'^rj =0 , 



(5) 



thus leading to an expression for r^' in terms of the dimen- 
sionless truncation radius of the corresponding spherical King 
model and a first estimate of the critical value of e 



-(0) _ 3 Ao _ 3 „ 
2 ao 2 



.(0): 



243 1/, 



-Ao 



(6) 



(7) 



The first exp ression can also be written as bf^ = 2/3 (see 
ISpitzeill987l) . We recall that Ao = Ao(r,r) and that r,r = r,r(^). 
Therefore, the right-hand side of Eqs. (|6]l and (|7]i depends only 
on the value of ^P. 

Here one important comment is in order Strictly speak- 
ing, the complete solution for a$c derived by the method of 
matched asymptotic expansions in § 4 of Paper I is a well- 
justified global uniform solution only for sub-critical (under- 
filled) models. Close to the condition of criticality, i.e., when 
r,r ^ fj, in the vicinity of the boundary surface the tidal term 

' We recall from Paper I that Aq = r,^^i/'Q°"'(rfr) and = Xg/frr, with 
^/iq"" the zeroth-order term of the asymptotic series of the internal solution. 
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Fig. 2. — Critical values of the extension parameter for first-order mod- 
els (dotted) and second-order models (solid) with different potentials of the 
hosting galaxy (u = 3,2, 1 from top to bottom). The dashed horizontal line 
shows the value 2/3 that is found when a zeroth-order approximation for a^c 
is used (see Eq. ^6)). 



eT (which is considered a small correction in the construc- 
tion of the asymptotic solution) becomes comparable to the 
cluster term a^c, so that the asymptotic solution is expected 
to break down. For such models, the iteration method de- 
scribed in § 5.2 of Paper I, which does not rely on the assump- 
tion that the tidal term is small, is preferred and expected to 
lead quickly to more accurate solutions. In practice, in line 
with previous work on the similar problem for rotating poly- 
tropes mentioned in the Introduction, we argue that the use 
of the second-order asymptotic solution constructed in Paper 
I will give sufficiently accurate solutions in the determination 
of the critical value of the tidal parameter e^,- from the sys- 
tem in Eq. (|4]i and in the consequent assessment of the gen- 
eral properties of models even when close to the critical case. 
The main reason at the basis of this argument is that, even for 
close-to-critical models, only very few stars populate the re- 
gion where the asymptotic analysis breaks down, so that the 
overall solution should be only little affected. A direct com- 
parison between selected critical models calculated with both 
the perturbation and the iteration method is presented in Ap- 
pendix A. 

If we make use of the full second-order asymptotic solution 
for the escape energy (^/;*"")*^^', in which the second-order ex- 
pressions for the cluster potential (recorded in Eq. ( IB3l l) and 
the cut-off constant are used, the system in Eq. (|4]l can be re- 
arranged and written in standard form 



with 



A(rr)e?. + 2B(rr)e., + C(rr) = 
D(rT)el- + 2E(rT)ecr+F(rT)-- 

3 [5 




--0 



(8) 



Airj) = Xifj- (bio - ^22^3) j\ -i-T + bAO 



aH, 



o = cto 



-b Al- 



ls 



+ bAA-7 
TT 16 



4 V TT 

15 /35 



45 



(9) 
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3 / 5 

Birr) = A 1 4 - (020 - 022^3) -\ -rl + 9iyrl , 

4 V TT 

CiPr) = 2Xofj , 



D(rr) = - + (^20 - /722\/3) 7 W - 4 

4 V TT 

^ 9 3 /S" 3 /35 

loy/n 8 V TT 16 V TT 



(10) 
(11) 



(12) 



(13) 



1/5 9 

£■(^7-) = air^ - Ai4 + (fl20-fl22%/3)- \/ -4 + ^ '^''r 

4 V TT / 

F(rr) = 2(ao4-Ao4), (14) 

where the relevant constants (A,- with / = 0,1,2, ai„ with 
I = 0,2 and bi,,, with / = 0,2,4 m = 0,2, ...,l), which are deter- 
mined by the matching process, are defined in § 4.4 of Paper 
I. [The corresponding system based on the first-order solution 
for V'*"" can be recovered by setting A{rx) = D(rj) = 0.] This 
system has been solved numerically, by means of the Newton- 
Raphson method, since the equations are nonlinear in Pt (in 
particular, they are polynomials of fifth and seventh-order, for 
the first and second-order solution respectively). As noted in 
the discussion of the simpler Eq. (|5]l, the solution for ecr can 
then be represented as a function of the concentration param- 
eter ^. The parameter space of the first and second-order 
models has been explored by means of an equally spaced grid 
from e = 5 X 10"^ to e = 1 x 10"^ at steps of 5 x 10"^ and from 
* = 0.1 to * = 10 at steps of 0.1. 

The parameter space for the second-order models is pre- 
sented in Fig. [T] The plot provides the contour levels of the 
extension parameter 5, with the uppermost solid line corre- 
sponding to S = 6cr ~ 2/3 (thus identifying the critical mod- 
els), based on the choice i' = 3 (Keplerian host galaxy). The 
critical curves for i' = 2 (host galaxy characterized by flat ro- 
tation curve) and for v = I (Plummer potential evaluated at 
Rq = b/ V2, with b the model scale radius) are shown as a 
dashed and dot-dashed line, respectively. 

Sub-critical, underfilled models (bottom-left corner of the 
figure), with 6 ^ (5„-, are only little affected by the tidal per- 
turbation. The maximally deformed models are those with 
5 « Scr (close to the uppermost solid line, i.e. close-to-critical 
configurations). Figure [T] shows that the critical value for the 
tidal strength parameter depends strongly on concentration, 
with a variation of almost four orders of magnitude in the ex- 
plored range of The figure also indicates that for lower 
values of v the critical curve moves upwards, i.e. the avail- 
able parameter space increases. 

The difference between the critical value of the tidal 
strength parameter for first and second-order models (for 
a chosen value of ly) is very small, around 10"^ for low- 
concentration models, down to 10"^ or less for models with 
5* w 10. The critical value of the extension parameter S de- 
pends only weakly on concentration and on as illustrated in 
Fig.E] 

In closing this section, se should reiterate that, in spite of the 
abundant use of symbols required by the analysis, the family 
of models that we have studied is characterized by two dimen- 
sionless parameters ('!';£). [As an alternative pair, we may 
refer to the standard concentration parameter C = log{r,r/ ro), 
equivalent to and frequently used in the context of spher- 
ical King models, and to the extension parameter 6 = r,r/rT, 
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Fig. 3. — Intrinsic density profiles (normalized to the central value) for 
critical second-order models with u = 3 and >!/ = 1 , 2, .... 10 (from left to right). 
Top panel (a): profile of the triaxial models along the i-axis (solid) and of the 
corresponding spherical King models (dashed). Bottom panel (b): profile of 
the triaxial models along the y-axis (solid) and the z-axis (dashed). 



equivalent to e.] The free constants A and a that appear in 
the distribution function fK{H) set the two physical scales. In 
turn, the models constructed by iHeggie & Ramamanil (Il995h 
are a one-parameter family of models, because these authors 
focused on the critical case and did not discuss the sub-critical 
regime. 

3. INTRINSIC AND PROJECTED DENSITY 
DISTRIBUTION 

3.1. Intrinsic density profile 

The models are characterized by reflection symmetry with 
respect to the three natural coordinate planes. With respect to 
the unperturbed configuration (i.e. the spherical King model 
with the same value of ^f), they exhibit an elongation along 
the i-axis (defined by the direction of the center of the host 
galaxy), a compression along the £-axis (the direction perpen- 
dicular to the orbit plane of the globular cluster), and only a 
very modest compression along the y-axis. 

Models with 5 < 0.2, regardless of the value of ^t, are 
practically indistinguishable from the corresponding spheri- 
cal King models; significant departures from spherical sym- 
metry occur for models with S « 0.4 or higher In Fig. |3]a 
we show the density profile along the i-axis for a selection of 
critical second-order models with = 3 in comparison with 
that of the corresponding spherical King models; note that 
for a model with "if = 2 the elongation is akeady significant 
at Logip/po) w -4, while for a model with ^' = 8 a simi- 
lar elongation is reached only at much lower density levels 
(Logip/po) « -7). The corresponding profiles along the y- 
axis and the z-axis are given in Fig.[3]b. 

For completeness, we checked the dependence of our den- 
sity profiles on the potential of the host galaxy. Consistent 
with the general trends suggested by Fig. |2] the elongation 
along the i-axis and the compression along the z-axis for the 
models with ly = 3 turn out to be slightly weaker than for the 
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Fig. 4. — Central values of the polai' (eo; dashed) and equatorial (r)o; solid) 
eccentricities of the isodensity surfaces of second-order models with u = 3, 
= 1,3,5,7 (from top to bottom) and e S [0, ecr('I')]. Horizontal lines mark 
the maximum value of eg ™d tjq reached by the critical models (i.e., for 



models with smaller values of h'. 

Since the dimensionless density distribution of a model, 
identified by (^,e), is given by p = p[ipir)], where ip is the 
dimensionless escape energy and p is a monotonically in- 
creasing function that vanishes for vanishing argument (see 
Eq. (11) in Paper I), there is a one-to-one correspondence 
between isodensity and isovelocity surfaces, the latter being 
defined by the condition ip^'"'\r) = S, where 5 is a constant 
(with < 5 < ^). We recall that nonspherical models of- 
ten exhibit equipotential surfaces rounder than the isodensity 
surfaces (e.g., see Evans 1993; Ciotti & Bertin 2005). The 
reason for the presence of this property in our models is that 
the supporting distribution function depends only on the Ja- 
cobi integral, i.e. the (isolating) energy integral in the rotating 
frame. Therefore, for each value of S, we can define the semi- 
axes a, b, and c of the corresponding triaxial isodensity sur- 
face, by means of the intersections of the surface with the i, 
y, and z axes, which turn out to follow the ordering a>b>c. 
The shape of the triaxial configuration can thus be described 
in terms of the polar and equatorial eccentricities, defined as 
e= [l-(c/a)2]'/^ and rj = [1 -(^/a)^]'/^, respectively. 

A surprising result can be derived analytically. In the inner- 
most region r <C Hr (i.e., for S ^ ^t), the dimensionless escape 
energy can be expanded to second order in the dimensionless 
radius 



s2 



--(l-zy)+A2oi'2o(e,</))+ (15) 



A22Y22(0,(j))] + -r^[(l+B2o)Y2o(e,^) + (l+B22)Y22(0,<j))] . 

Here some terms of the second-order solution do not con- 
tribute (e.g., it can be readily checked that V"! 4m('^) ^ ''^ ™^ 
ip^^ir) - r*-). Then by setting V'^'"'\fl,0,0) = V'^'""(0,^,0) = 
V'*""^(0,0,c), we find that in the innermost region the eccen- 
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Fig. 5. — Profiles of the polar (e; dashed) and equatorial (rj; solid) eccen- 
tricities of the isodensity surfaces for selected critical second-order models 
with u = 'i and >!/ = 1 , 3, 5, 7 (from left to rig ht). Dotte d horizontal lines show 
the central eccentricity values (see Eqs. )161 and )17t ). 



tricities tend to the following non-vanishing central values 

{ e(A22 - \/3A2o + (e/2) [ 1 + B22- V3( 1 + 820)])^/^}''^ 



{6 + 2e[3( 1 - i.) - A20 - 1 + B20) '''' 

(16) 



{e[2A22 + (l+g22)£]v/T5A}'/' 
{6 + et/i + (eV2)£/2}'/2 



where 



t/i = 6( 1 - 1/) + (A20 + \/3A22) ^^5/^ 
d2 = [1 +B20+ \/3(l +B22)]j5j^ 



(17) 



(18) 
(19) 



which depend explicitly on the tidal strength and implicitly 
on the concentration. This result is nontrivial. In fact, since 
the tidal potential is a homogeneous function of the spatial 
coordinates, naively we might expect that in their central re- 
gion the models reduce to a perfectly spherical shape (i.e., 
eo = '70 = 0), even for finite values of the tidal strength. In- 
stead, eo and 770 are ©(e'/^) and strictly vanish only in the 
limit of vanishing tidal strength. 

Figure |4] shows the central values of the eccentricities for 
second-order models with v = !> and selected values of con- 
centration, as a function of tidal strength within the range 
[0, ecri"^)]. Consistent with the general trends identified in the 
discussion of the parameter space, low-concentration models 
show the most significant departures from spherical symme- 
try. The full eccentricity profiles (as a function of the major 
axis) are shown in Fig. |5] for a selection of critical second- 
order models; here the calculation of e and 77 has been per- 
formed by numerically determining the values of the semi- 
axes of a number of reference isovelocity surfaces, defined 
by i^'"'\f) = Si = (25-0^'/25-0.01 with / = 0, ..,25. Outside 
the central region, the profiles increase monotonically and, in- 
dependently of concentration, at the boundary they reach ap- 
proximately a fixed maximum value (e w 0.78 and rj « 0.74), 
which corresponds to the fact that the shape of the boundary 
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surface of a critical model ((5J;°^ = 2/2>; see also Fig.|2]i depends 
only modestly on concentration. 



3.2. Compa rison with the models constr ucted by 
\Hessie & Ramamani ( 179931) 

The method used in Paper I for the construction of the mod- 
els illustrated in this paper can be summarized as follows. The 
solution in the internal (Poisson) and external (Laplace) do- 
mains are expressed as an asymptotic series with respect to 
the dimensionless parameter e, representing the tidal strength 
(defined in Eq. (|2|i), which is considered to be small. The 
A:th-order term of the asymptotic series of the internal and ex- 
ternal solution are denoted by ip'l.""(r) and respectively, 
so that the zeroth-order terms define the standard spherical 
King models (see § 4.1 in Paper I for a detailed discussion). 
The quantity (^^'"'>)<*) indicates the A:th-order external solution, 
i.e. the corresponding asymptotic series truncated at the term 
ipf''^; a similar notation holds for the internal solution. The 
validity of the expansion breaks down where the second term 
is comparable to the first, i.e. where ^jjo = 0(e). This sin- 
gularity is cured by introducing a boundary layer in which 
both the spatial coordinates and the solution are suit- 
ably rescaled with respect to the tidal parameter To obtain a 
uniformly valid solution over the entire space, an asymptotic 
matching (see Van Dvke 1975; Eq. (5.24)) is performed be- 
tween the pairs (?/;<'"",'(/'""=") and (-^'''''W'"")- Each term ipk(f) 
is then expanded in spherical harmonics with radial coeffi- 
cients ipkjmir). The internal region requires a numerical so- 
lution of the Cauchy problems for the radial coefficients (we 
used a fourth-order Runge-Kutta code) while in the external 
region a formal solution with multipolar structure is available 
and in the boundary layer the integration in the radial variable 
can be performed analytically (see § 4.2 and 4.3 in Paper I, 
respectively). 

The models described by iHeggie & Ramamanil (Il995h are 
also based on a perturbation approach, but the method used 
is different from ours, provides a solution of the Poisson 
equation that is first-order with respect to the tidal parame- 
ter, and is restricted to the "critical" case; actually, as noted 
in § 2 after Eq. (|7]i, the perturbation approach is bound to 
break down in the critical case. Technically speaking, their 
method is in the form of a "patching" procedure, in contrast 
with our asymptotic matching. Therefore, the models con- 
structed in Paper I, while consistent, to first order, with those 
of Heggie & Ramamani (1995), are more general. We also 
recall that our method is also applicable to systems described 
by different distribution functions (see Appendix B and C in 
Paper I). 

We have thus performed a quantitative comparison between 
the intrinsic density profiles of our c ritical second-order mod - 
els and those of the models- by Heggie & Ramamani (1995J), 
both referred to the case in which the host galaxy is Keplerian 
(v = 3). As desired, there is substantial consistency except for 
the outermost part of the boundary layer in which our models 
are slightly more compact, due to a global "boxiness" effect 
induced by the second-order term present in our models in 
which harmonics of order I = 4 also play a role. In Fig. |6]we 
represent the relative difference between the two density pro- 

^ For this purpose, we used t he code that implements the models described 
bv lHeggie & Ramamanil II 9951) . writt en by D. C. Heggie and availa ble within 
the STARLAB software environment jPortegies Zwart et alJlOOlT) 
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Fig. 6. — Relative difference between the intrinsic density profiles of 
critical second-order models (p;;) constructed in this pap er and those of the 
corresponding (first-order) models (phr) described by Heggie & Ramaman| 
1 1995). The comparison has been performed along the three axes in the whole 
intemal+boundary region (see main text). Here we illustrate the difference 
along the ,r-axis in the boundary layer ('I' = 1,2,..., 10 from left to right). 
At variance with Fig. |3] the spatial coordinate is scaled with respect to the 
truncation radius instead of the scale radius tq. 

files evaluated along the x-axis (with the coordinates scaled 
with respect to the truncation radius instead of the usual scale 
radius ro) for selected values of "if. A similar behaviour is 
found also along the y-axis, for y/r,,- in the range [0.80,1], 
and along the z-axis, for z/r,,- in the range [0.80,0.95]. In the 
central part of the internal region, along the principal axes, the 
relative difference is smaller than 5 percent for every value of 
5* we tested, while near the transition to the boundary layer 
(i.e x/r,r < 1 and y/r,r, z/r,r < 0.8) a difference of 20 per- 
cent can be reached in the case of highly concentrated mod- 
els (^f = 8,9, 10). We interpret these differences as due to 
the combined effects of the patching vs. matching adopted 
process and of the different grid on which the Cauchy prob- 
lems for the radial coefficients are solved (we used a regu- 
lar radial grid while Heggie & Ramamani ( 1995) used a more 
complex tabulation resulting from their choice of taking the 
zeroth-order cluster potential as the independent variable and 
of the function /n(l + r^) instead of r). 

3.3. Global quantities 

The previous discussion has focused on the shape of the iso- 
density surfaces of the models. In particular, some interesting 
conclusions have been derived based on a local analysis of 
the central region and of the outer boundary of the configu- 
ration. The maximal departures from spherical symmetry are 
reached at the periphery, but these hinge on the distribution 
of the very small number of stars that populate the outer re- 
gion of the cluster We may thus wish to study some global 
quantities that better characterize whether significant amounts 
of mass (and, correspondingly, of light) are actually involved 
in the deviation of the model from spherical symmetry. One 
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Fig . 7 . — Ratio of two pairs of quadrupole moments for the second-order models with ^ = 5,0.1 < <5 < ScA'^), and v = 1,2,3 (in the left panel, the sequence 
of models with different v values runs from bottom to top; in the right panel, it runs from top to bottom). Th e val ues obtained from numerical integration over 
the entire triaxial volume (dots) are compared to the analytical approximations (solid line) given by Eqs. (26) and (25); the analytical estimates of the ratios for 
first order models are also shown (dashed horizontal lines). The propagation of the errors of the numerical integration leads to the plotted error bars. 



Standard such global measure is provided by the quadrupole 
moment tensor 

Qij= I (3xiXj-r^Sij)p(r)d^r = 
Jv 



An 



/ (3i(ij- 

V 



-f'5ij)p{r)d^r = ArlQ, 



(20) 



with the integration to be performed in the volume V of the 
entire configuration. Here the notation for the function p and 
for the constant A is the same as in Eq. (11) of Paper I. In con- 
trast with the frequently used inertia tensor I/j = JypxiXjd^r 
(e.g., see Chandrasekhar 1969; chap. 2), the quadrupole mo- 
ment is defined in such a way that in the spherical limit it 
vanishes identically. In our coordinate system it is diagonal. 
Note that the tidal distortions require that the non-vanishing 
terms of the inertia tensor follow the ordering 7^^ > A v ^ hz'^ 
to visualize the geometry of the system, we may thus also re- 
fer to the average polar and equatorial eccentricities e and fj 
defined by the relations lyy = (1 - fj^)Ixx, Az = (1 - e^)Ixx- In 
general, we have Qyy/Qxx = {e^-2rf)/{e^ + ff-) and Qzz/Qxx = 
{ff - 2e^) /{e^ + ff ), with the prolate configuration identified 

by e = fy, i.e. Qyy/Qxx = Qzz/Qxx = "1/2. 

Since most of the mass is contained in the inner regions, 
global quantities can be evaluated approximately by neglect- 
ing the contribution from the region corresponding to the 
boundary layer. We can thus use the second-order solutions 
for p obtained in Paper I by the method of matched asymptotic 
expansions and conveniently reduce the calculation of global 



quantities to an easier integration in spherical coordinates in- 
side the sphere of radius r,r. Therefore, for the quadrupole 
moment tensor we find 



,(2) . 



G/>,ie + 2/^.2 Y 



(21) 



We emphasize that this estimate is expected to be a good ap- 
proximation only for those models for which the contribution 
of the boundary layer is negligible with respect to the one of 
the internal sphere of radius r,r- 

The relevant components on the diagonal can be expressed 
in terms of the matching constants of the external solution (for 
the relevant definitions, see Eqs. (63) and (69) in Paper I) 

A" 



-VSnpi^) 



aioe + bio— 



V3 ( 0226 + ^22- 



^2\ 1 



(22) 



a2oe + b20— 



Q 



i(2) . 



-VSwpi^) 

+ V3 ^fl22e + ^22Y^ 



4 ^ / e2 
--VSnpC^) [a2Qe + b2o— 



(23) 



(24) 



We recall that the constants 020 and /720 are positive, while 
fl22 and ^722 are negative (and larger in magnitude). Therefore, 
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Q^^^ is positive and 2*^? and Ql? are negative, consistent with 
the detailed elongation and compressions observed in the den- 
sity profile. A summary of the derivation of these formulae is 
provided in Appendix B. 

As a measure of the degree of triaxiality of a given config- 
uration, we have calculated the following ratios 



g(2) 



g: 



yy 

,(2) 



(a2() + ^2oe/2) + \/3 (fl22 + ^22e/2) 



(a20 + ^2oe/2)-\/3 (fl22 + ^22e/2) 

J,, _ -2 (fl20 + ^2oe/2) 

W~ (fl20 + W2)-\/3(fl22 + W2) ' 



(25) 



(26) 



which depend explicitly on the tidal strength parameter and 
implicitly on the concentration parameter. In the limit of van- 
ishing tidal strength, we find 



Go 
Go 



_ T20(rn) + V3T22(rn) 



1 



,(1) 



G 



2v+\ 
2 + v 



Gix T2Q(f,r)-V^T22(f,r) 2l^+l 



(27) 



(28) 



where T2m(r) are the quadrupole coefficients of the tidal po- 
tential (see Eqs. (38) and (39) in Paper I). This result is non- 
trivial, because, in this limit, numerator and denominator are 
both expected to vanish. Note that only for h' = I the ratio 

Gvv/Gxv = 0(e). 

Earlier in this paper we mentioned that two physical scales, 
such as total mass and central velocity dispersion, correspond 
to the two dimensional constants A and a that appear in the 
distribution function fxiH). In fact, the total mass of the sys- 
tem is given by 



M= / p{r)(Pr = Arl 



p(r)d^r = Ar^M 



(29) 



If we insert the second-order solution for p obtained in Paper 
I, we find 







/ 


b [ dOsine [ 


Jo 


Jo Jo 



?2 



with 



: Mq + eMi + y -^2 



4Trp(^) 
Mi = Mi(^) = ^— ^A/ 



(30) 



(31) 



Here each term of the expansion is related to the correspond- 
ing constant with I = (i.e., the monopole term) of the expres- 
sion of the external solution of the Poisson-Laplace equation 
calculated by means of the method of matched asymptotic ex- 
pansion (for the relevant definitions, see Eqs. (59), (61), (67) in 
Paper I). 

The quality of the analytical estimates for the total mass and 
the quadrupole moment tensor has been checked by compar- 
ing the values obtained from asymptotic analysis with those 
resulting from direct numerical integration of Eq. ( |29] l and 
Eq. (l20b respectively, in which the density profile p = /5[?/;(r)] 
is used without any additional expansion. The integration of 
the distribution function over the entire space, required by 
those global quantities, has been performed by means of an 
Adaptive M onte Carlo method (the algorithm VEGAS, see 
iPress et an[l992 : § 7.8), well suited for our geometry. For the 
quadrupole, the results are illustrated in Fig.|2l For the mass. 
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Fig. 8. — Average polar (e: dashed) and equatorial (77; solid) eccentricities 
for critical second-order models as a function of concentration, with 0. 1 < 
< 10, for two different host potentials. Models with u = 3 con'espond 
to the inner pair of curves, while models with u = i to the outer pair. The 
eccentricities have been determined numerically from their definitions (see 
main text) in terms of p = /3[t/>(f)], without any additional expansion; the 
relative errors are around 10"^. 

the dimensionless function M*^'(vl') is basically unchanged 
(within 0.5 percent) with respect to the function character- 
izing the spherical King models; the Monte Carlo integration 
is very accurate, with relative errors around lO"'', and the an- 
alytical approximation given by Eq. (|30] | shows an excellent 
agreement for every value of in the whole range of the ex- 
tension parameter [0, Scri'i')]- 

The average eccentricities for critical second-order mod- 
els as a function of concentration, with 0.1 < ^ < 10, for 
two different choices of the host galaxy potential (1/ = 3,1), 
are shown in Fig. [8] The values are calculated directly from 
the definitions given earlier in this subsection in terms of p = 
p[ilj(r)] with no additional expansions. A non-monotonic de- 
pendence on concentration is revealed, with generally higher 
average eccentricities attained by low-concentration models. 
The trends of the polar and equatorial eccentricities are basi- 
cally the same, as shown by the fact that the related curves 
in the plot can be matched approximately by a rigid trans- 
lation. As expected (see § 3.1), models with u = I show a 
larger separation between polar and equatorial eccentricities 
than models with v = 3. The presence of a minimum for the 
curves at « 6.5 occurs approximately at the location where 
the function Log[ecri'^)] shows an inflection point (regardless 
of the value of v; see Fig.lTJ. 

3.4. Projected density profile 

Under the assumption of a constant mass-to-light ratio, pro- 
jected models can be compared with the observations. We 
have then computed surface (projected) density profiles and 
projected isophotes. The projection has been performed along 
selected directions, identified by the viewing angle [9, </)), cor- 
responding to the zp axis of a new coordinate system related to 
the intrinsic system by the transformation i/) = Rx; the rotation 
matrix R = Ri{d)RT,{(f)) is expressed in terms of the v iewing 
angles, by taking the xp axis as the line of nodes (see iRvdenI 
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Fig. 9. — Projections of a second-order critical model = 2 and u = 3) along the x axis (left panel [a]) and along the y axis (central panel [b]). The ellipticity 

profiles (light panel [c], from bottom to top) refer to the projections along directions identified by (6 = n/2, <f) = in 16) with i = 3; dots represent the locations 

of the isophotes drawn in panels [a],[b], which correspond to selected values of S/So in the range [0.9, 10"*]. The an'ow indicates the position of the half-light 
isophote (practically the same for every projection considered in the figure). 



119911; for an equivalent projection rule). Given the symmetry 
of our models (see § 3.1), viewing angles can be chosen from 
the first octant only. In particular, we used a 4 x 4 polar grid 
defined by 6',- = / 7r/6 and (pj = j 7r/6 with /, j = 0, . . , 3, and cal- 
culated (by numerical integration, using the Simpson rule) the 
dimensionless projected density 



S](ip,yp) = / p(rp)dzp , 



(32) 



where Zsp = (x^-Xp-ypy^^ with the edge of the cluster 
along the i axis of the intrinsic coordinate system (i.e., we 
"embedded" the triaxial configuration in a sphere of radius 
given by its maximal extension). The projection plane (xp,yp) 
has been sampled on an equally-spaced 108 x 108 square grid 
centered at the origin (note that for Xp+yj,> x^ the projected 
density is correctly set to zero). 

The morphology of the isophotes of a given projected im- 
age can be described in terms of the ellipticity profile, defined 
as e = 1-bp/ap where ap and bp are the semi-axes, as a func- 
tion of the major axis dp. As already noted for the (intrinsic) 
eccentricity profiles, the deviation from circularity increases 
with the distance from the origin. In the inner region, the el- 
lipticity is consistent with the central eccentricities eo and ryo 
calculated in § 3.2. 

By taking lines of sight different from the axes of the sym- 
metry planes, we have also checked whether the projected 
models would exhibit isophotal twisting. For all the cases 
considered, the position angle of the major axis remains un- 
changed over the entire projected image. Tests made by 
changing the resolution of the grid confirm that, even in the 
most triaxial case (v = 1), no significant twisting is present. 

The first two panels of Fig. |9] show the projected images 
of a critical second-order model with ^ = 2 and v = 3 along 
the (7r/2,0) and (7r/2,7r/2) directions, (i.e., the x and y axis of 
the intrinsic system), corresponding, respectively, to the least 
and to the most favorable line of sight for the detection of 
the intrinsic flattening of the model. For the same model, the 
third panel illustrates the ellipticity profiles for various lines 
of sight. 

Figure [TO] shows the surface density profiles along the xp 
and yp axes of the projection plane for ten critical second- 
order models with i' = 3, viewed along the {tt/2,tt/2) direc- 
tion. As a further characterization, for the same models in the 



lower panel we also present the surface density profiles ob- 
tained by averaging the projected density distribution on cir- 
cular annuli; this conforms to the procedure often adopted by 
observers in dealing with density distributions with very small 
departures from circular symmetry (e.g., see iLanzoni et al.l 
|2007). As expected, circular-averaged profiles lie between the 
corresponding regular profiles taken along the principal axes 
of the projected image. 

4. INTRINSIC AND PROJECTED KINEMATICS 

By construction, the models are characterized by isotropic 
velocity dispersion. The intrinsic velocity dispersion can be 
determined directly as the second moment in velocity space 
(normalized to the intrinsic density) of the distribution func- 
tion 



.^W=^4M = V(,), 



5fl7(5/2,7^) 



(33) 



where 7 represents the incomplete gamma function (near the 
boundary of the configuration, the velocity dispersion profile 
scales as (T^it/j) ^ {2/l)ip). This shows that the isodensity sur- 
faces of the models are in a one-to-one correspondence with 
the isovelocity and isobaric surfaces (defined by <J^[il^(r)] = 
const). As noted for the intrinsic density profiles in § 3.1, 
a compression along the z axis and an elongation along x 
axis occur also for the intrinsic velocity dispersion profiles. 
In Fig. [TT]a we present the intrinsic velocity dispersion pro- 
files along the i axis for the same critical models illustrated in 
Fig.[3]compared to the profiles of the corresponding spherical 
King models. The behavior of the projected velocity disper- 
sion profiles near the boundary is significantly different from 
that of the spherical models. 

The projected velocity moments can be calculated by in- 
tegrating along the line of sight (weighted by the intrinsic 
density) the corresponding intrinsic quantities. Therefore, the 
projected velocity dispersion is given by 



(Tp(xp,yp)-- 



Xfr aHrp)p(rp)dzp 



^(xp,yp) 

2 X!r 7[7/2,V(fp)]exp[^(fp)]affp 



5a 



Y.(xp,yp) 



-ajixpjp) . (34) 
a 
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Fig. 10. — Projected density profiles (normalized to the central value) for 
the same ten second-order critical models displayed in Fig. [3] Top panel (a): 
the models are viewed from the y axis and the profiles taken along the two 
principal axes in the projection plane (along icp (solid) and yp (dashed), i.e. 
along the ,v and z axes of the intrinsic frame of reference). Bottom panel (b): 
the projection is performed on the same line of sight of the previous panel, but 
the profiles are taken by averaging the projected surface density on circular 
annuli, as if the image were intrinsically circular 

Figure [TT]b shows the projected velocity dispersion profiles 
along the xp and yp axis of the projection plane for the same 
models displayed in Fig. [10] (the line of sight is defined by 

(7r/2,7r/2)). 

5. DISCUSSIONS AND CONCLUSIONS 

5.1. A comment on the complex structure of the outer 
regions 

In Paper I we have noted the singular character of the math- 
ematical problem associated with the free boundary set by the 
three-dimensional surface of these tidally truncated models. 
In § 2 of the present Paper, we have further emphasized the 
additional singularity that characterizes close-to-critical mod- 
els (see comment after Eq. (Q, which prompted the analysis 
described in the Appendix A). As is true in general in the 
study of boundary layers and similar problems, it is no won- 
der that in the vicinity of such critical boundaries, a number of 
complex physical effects may take place and play an impor- 
tant role in determining the detailed structure of the solution. 
On the other hand, the properties of the derived solution away 
from the boundary are quite robust (see § 3.3 and the Ap- 
pendix). As to some of the subtle properties of the expected 
distribution function close to the edge of a cluster, it is inter- 
esting to summarize here below the main results that emerge 
from a vast body of evolutionary models (N-body, Fokker- 
Planck, Monte Carlo, gas), on the issue of the interplay be- 
tween pressure anisotropy and tidal effects. 

Since the first solutions of the Fokker-Planck equa tion by 
means of a Monte Carlo approa ch, as described by iHenonI 
(1197 lb or ISpitzer & Harj (119711) . it has been shown that 
one-component isolated spherica l clusters, sta rting from a 
variety of initial conditions (see ISpitzerl [1987^, chap. 4 and 
references therein), develop a core-halo structure in which 
the central parts are almost isotropic while the outer re- 



gions are characterized by radial anisotropy. A com- 
monly reported interpretation is that the halo is mainly 
populated by stars scattered out from the core on radial 
orbits. If the evolution of a cluster is influenced not 
only by internal processes but also by the presence of 
an external tidal field, the growth of pressure anisotropy 
is signifi cantly damped. Direct N-body simu l ations 
(e.g., see iG iersz & Heggiel 119971; l Aarseth & Heggiel 119981; 
iBaumgardt & Makino 20 03|: [Leee t al. 2006), anisotropic 
Fokker-Planck (e.g., see Takahashi et al.lll997l) . and Monte 
Carlo models (e.g., see ,Gier sz 2001.) of both single and multi- 
mass systems suggest t hat clusters in circular orbits (and even 
in eccentric orbits, see IBaum gardt & Makinoll2003l) . starting 
from isotropic initial conditions, tend to preserve pressure 
isotropy, except for the outermost parts which become tan- 
gentially anisotropic due to the preferential loss of stars on ra- 
dial orbits, induced by the tidal field. The overall agreement 
on this result is nontrivial, because of the aforementioned dif- 
ferences in the treatment of the external tidal field. Even 
extreme cases of time-dependent tides, such as disk shock- 
ing, influence the degree of pressure anisotropy since it has 
been shown that they may represent a dominant mechanism 
("shock relaxation") of the energy redistribution, le ading to a 
substantial isotropy, of the weak ly bound stars (see lOh & LinI 
[1991 iKundic & OstrikeiifT995i ; both papers are based on a 
Fokker-Planck approach). 

These results confirm that, of course, the properties of the 
models constructed in Paper I and described in this paper 
should be taken only as a zeroth-order reference frame, to 
single out the precise role of external tides, and should not 
be taken literally as a realistic representation of real objects 
since a number of simplifying assumptions are made. On the 
other hand, by comparing data with such an idealized refer- 
ence model it will be possible to better assess the role of tides 
with respect to other physical ingredients studied separately. 

5.2. Summary and concluding remarks 

The main results that we have obtained from the detailed 
analysis of the family of tidal triaxial models can be summa- 
rized as follows: 

1. Two tidal regimes exist, namely of low and high- 
deformation, which are determined by the combined 
effect of the tidal strength of the field and of the 
concentration of the cluster The degree of deforma- 
tion increases with the degree of filling of the rele- 
vant Roche volume. Far from the condition of Roche 
volume filling, the models are almost indistinguish- 
able from the corresponding spherical King models. A 
number of studies have investigated the evolution of 
tidally perturbed stellar s ystems initially underfilling 
their Roche lobe (e.g ., see iGieles & Baumgardtil200"8c 
IVesperini et al.ll20()9l) . concluding that some of the rel- 
evant dynamical processes, in particular evaporation, 
depend on the degree of filling of the Roche volume. 
The intrinsic properties of the models discussed in this 
paper can be useful for setting self-consistent nonspher- 
ical initial conditions of numerical simulations aiming 
at studying in further details the evolution of configura- 
tions starting from sub-critical tidal equilibria. 

2. For a given tidal strength, there exists a maximum value 
of concentration for which a closed configuration is al- 
lowed (see Fig. [Tjl. In such "critical" case, the trunca- 
tion radius of the corresponding spherical King models 
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Fig. 1 1 . — Top panel (a): intiinsic velocity dispersion profiles (solid, nor- 
malized to the central value) along the f axis for the same ten second-order 
critical models illustrated in Fig.[3]a compared to the corresponding spheri- 
cal King models (dashed). Bottom panel (b): projected velocity dispersion 
profiles (normalized to the central value) for the same ten second-order crit- 
ical models illustrated in Fig. |10| a, viewed along the same direction. Solid 
and dashed lines show the profiles along the xp axis and yp of the projection 
plane, respectively. 

is of the same order of the tidal radius of the triaxial 
model. The shape of the boundary of the maximally de- 
formed models is given by the "critical" zero-velocity 
surface of the relevant Jacobi integral and is basically 
independent of concentration, while the deformation of 
the internal region strongly depends on the value of the 
critical tidal strength and is more significant for low- 
concentration models. This statement agrees with a 
general trend noted by White & Shawll ( 119871) for the 
globular cluster system of our Galaxy. 

3. The structure of the models can be described in terms 
of the polar and equatorial eccentricity profiles of the 
intrinsic isodensity surfaces. The maximal departures 
from spherical symmetry are reached in the outskirts. 

4. For finite tidal strength, the central values of the polar 
and equatorial eccentricities are finite, C9(e'''^); this re- 
sult is nontrivial since the tidal potential which induces 
the perturbation is a quadratic homogeneous function 
of the spatial coordinates. 

5. Global measures of the degree of triaxiality in terms 
of the quadrupole moment tensor have been introduced 
and calculated for different values of the tidal strength 
and different potentials of the host galaxy. The potential 
of the host galaxy sets the geometry of the tidal pertur- 
bation, as nicely shown by the analytic expressions for 
the ratios of the components of the quadrupole moment 



tensor, which, in the case of first-order models, reduce 
to simple functions of the coefficient v. Average ec- 
centricities have been calculated from the inertia tensor 
components, evaluated numerically. 

6. As a key step in the direction of a comparison with ob- 
servations, projected density profiles and ellipticity pro- 
files have been calculated for a number of models for 
several lines of sight. 

7. The study of the relevant (projected) isophotes indi- 
cates that no isophotal twisting occurs. This result 
is nontrivial since the models are nonstratified and 
centrally-concentrated, conditions under which, in prin- 
ciple, isophotal twist may occur (see IStarkI (Il977l) : 
models based on Stackel potentials are also known to 
be twist-free, as shown by lFranxl(ll988h ). 

8. Finally, close to the boundary, the intrinsic and pro- 
jected kinematics shows significant differences with re- 
spect to that of spherical models. 

Since our models are all characterized by monotonically in- 
creasing ellipticity profiles, they cannot explain the variety 
of be havior of observed ellipticity profiles (see iGeyer et al.l 
fT983h . but the range of the predicted flattening {e < 0.3) is 
consistent with that ob served in most globular clusters (see 
^ite & Shawl" (1987) for the clusters in the Milky Way and 
Frenk & Fall (1982) for those in the LMC). Therefore, with 
the exception of special clusters such as Omega Centauri, we 
think that the modest but frequently observed deviations from 
spherical symmetry might have their origin traced to tides. 

Finally, since our models are intrinsically more elongated 
(in the direction of the center of the host galaxy) than spher- 
ical King models, they might be useful for interpreting clus- 
ters with a surface brightness profile extending beyond that 
predicted by the spherical King models. Rec ent investiga- 
tions (see McLaughlin & van der Marell 12005 ^ suggest that 
such "extra-tidal" structures are a fairly generic feature, espe- 
cially for extragalactic clusters, and not just a transient prop- 
erty, present only at young ages. 

On the basis of the work presented here, in a following pa- 
per we plan to address in detail the issues raised by observa- 
tions. 
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APPENDIX 
A. PERTURBATION VS. ITERATION 

For completeness, we calculated selected models also by means of the iteration method described in § 5.2 of Paper I, in order 
to verify the quality of the solution obtained with the method of matched asymptotic expansions, in particular in the critical 
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regime. This technique follows the approach proposed by iPrendergast & Tomeii ( Il970h . The basic idea is to get an improved 
solution of the Poisson equation (see Eq. (77) in Paper 1) by evaluating the right-hand side with the solution obtained from the 
immediately previous step, starting from the spherical King models taken as initial "seed solutions". In our code, the required 
spherical harmonic analy sis and synthesis of density and poten tial have been perform ed by means of the S2kit 1.0 package 
(iKostelec & Rockmord l2004). which makes use of FFTW 3.2.1 (Frigo & Johnson 2005). We decided to truncate the harmonic 
series at / = 4 in order to be consistent with the maximum harmonic index admitted by the second-order asymptotic solution. The 
iteration stops when convergence to four significant digits in the whole domain of the solution is reached. 

For selected values of the concentration parameter in the range 0.1 < 5* < 10 (for simplicity, we considered only the case of an 
external potential with v = 3), the corresponding critical value of the tidal strength parameter obtained with the iteration method 
is consistent to 10"^ with the value determined by the numerical solution of Eq. (O, in which the constants obtained from the 
asymptotic matching are used. For a critical model up to 20 iteration steps are required for convergence, while a subcritical 
solution typically takes in 4 to 8 steps. 



B. GLOBAL QUANTITIES FROM THE MULTIPOLE EXPANSION OF THE CLUSTER POTENTIAL 

Based on the expansion in spherical harmonics of l/|r-r'| given in Eq. (3.70) of iJacksonI (Il999l) . the external potential 
generated by our model can be expressed by means of the following multipole expansion 

„ +00 , +/ 



a$5"'(r) = - 



P(*) 



^2/+l ^ 

/=0 m=-l 



lilt 



£1+1 



(Bl) 



47r/5(^') JV 

where the multipole coefficients are defined as 

q,„,= f d'r'Y,„,(0',c^')r"p(r'), 
Jv 

with the integration to be performed in the volume V with boundary surface defined hy = 0. This general expression can be 
compared with the second-order solution of the Laplace equation obtained in Paper I 



(B2) 



Ao + Aie+A2y 
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Y20(e,<l>) 



a22f- + b22 — 



Y22(0,<l>) 



740(^,0) , ^ e^ Y42(e,(f,) , ^ £^ 744(^,0) 

+t>40 TT + ''42 TT + ^^44 TT 



(B3) 



(2) 

in order to determine the relation between the second-order multipole coeffici ents q)J (i.e. calculated by means of the second- 

(B4) 



order expression for the density) and the matching constants that appear in Eq. ( IB3b . Therefore, we find that 

.2 9 



Ao + Aie+A2y =- 



-(2) 



/47r/5(*) 

so that Eq. (IJTT i follows. For the higher-order harmonics we find the following relations for the non-vanishing coefficients 

.2 9 



a2,ne + b2,„-= ^^^^ 



<2) 
l2m ' 



for m = 0,2 and 



form = 0,2,4. 

Recalling that we are using real spherical harmonics with the Condon- Shortley phase, we get 



1 



^20 = 



Qzz 



q22 = -rT\ (Gxx-fivv) 

12 V TT 
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(B6) 



(B7) 



(B8) 



To determine all the nontrivial components of the quadrupole moment tensor, we use the condition Tr(2,j) = 0. Therefore, for the 
second-order solution of Paper I, the system 



20 / TT 

20 



ei?-Gg^=-yA/T7/5W ( «22e + ^22- 



(B9) 



A(2) A(2) A(2)^Q 
V. WJxt ^yy ^zz " 

leads to the expressions recorded in Eqs. (I22]l- (l24b . 
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